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In  this  paper,  spatially  high-order  Residual-Distribution  (RD)  schemes  using  the  first- 
order  hyperbolic  system  method  are  proposed  for  general  time-dependent  advection-dif- 
fusion  problems.  The  corresponding  second-order  time-dependent  hyperbolic  advection- 
diffusion  scheme  was  first  introduced  in  [NASA/TM-2014-218175,  2014],  where  rapid  con¬ 
vergences  over  each  physical  time  step,  with  typically  less  than  five  Newton  iterations,  were 
shown.  In  that  method,  the  time-dependent  hyperbolic  advection-diffusion  system  (linear 
and  nonlinear)  was  discretized  by  the  second-order  upwind  RD  scheme  in  a  unified  manner, 
and  the  system  of  implicit-residual-equations  was  solved  efficiently  by  Newton’s  method 
over  every  physical  time  step.  In  this  paper,  two  techniques  for  the  source  term  discretiza¬ 
tion  are  proposed;  1)  reformulation  of  the  source  terms  with  their  divergence  forms,  and 
2)  correction  to  the  trapezoidal  rule  for  the  source  term  discretization.  Third-,  fourth,  and 
sixth-order  RD  schemes  are  then  proposed  with  the  above  techniques  that,  relative  to  the 
second-order  RD  scheme,  only  cost  the  evaluation  of  either  the  first  derivative  or  both  the 
first  and  the  second  derivatives  of  the  source  terms.  A  special  fourth-order  RD  scheme  is 
also  proposed  that  is  even  less  computationally  expensive  than  the  third-order  RD  schemes. 

The  second-order  Jacobian  formulation  was  used  for  all  the  proposed  high-order  schemes. 

The  numerical  results  are  then  presented  for  both  steady  and  time-dependent  linear  and 
nonlinear  advection-diffusion  problems.  It  is  shown  that  these  newly  developed  high-order 
RD  schemes  are  remarkably  efficient  and  capable  of  producing  the  solutions  and  the  gra¬ 
dients  to  the  same  order  of  accuracy  of  the  proposed  RD  schemes  with  rapid  convergence 
over  each  physical  time  step,  typically  less  than  ten  Newton  iterations. 

I.  Introduction 

In  this  paper,  we  construct  very  efficient  high-order  schemes  for  general  time-dependent  advection- 
diffusion  problems,  based  on  the  first-order  hyperbolic  system  method.  ’  In  this  method,  the  diffusion  term 
is  reformulated  as  a  hyperbolic  system,  leading  to  the  unification  of  advection  and  diffusion  as  a  single  hyper¬ 
bolic  system.  The  drastic  change  in  the  type  of  equations,  parabolic  to  hyperbolic,  brings  several  dramatic 
improvements  in  the  construction  of  numerical  schemes:  hyperbolic  schemes  for  diffusion,  the  same  order  of 
accuracy  for  the  solution  and  the  gradients,  orders-of-magnitude  convergence  acceleration,  etc.,  which  have 
been  demonstrated  for  steady  diffusion  and  viscous  problems  in  Refs.  “  and  unsteady  advection-diffusion 
problems  in  Ref.  1  It  is  based  on  the  reformulation  of  the  governing  equations,  and  therefore  applicable  to 
any  discretization  method.  In  this  work,  we  consider  a  Residual-Distribution  (RD)  method,  which  has  been 
well  developed  for  hyperbolic  systems  and  has  a  superior  feature  of  achieving  second-order  accuracy  in  a 
compact  stencil. 

In  the  previous  work,  we  extended  the  hyperbolic  method,  for  the  first  time,  to  time-accurate  com¬ 
putations  by  an  implicit  time-integration  method  based  on  the  second-order  backward  difference  formula. 

‘Research  Aerospace  Engineer,  Aerothermodynamics  Branch,  M/S  408A,  Senior  AIAA  Member,  Ali.R.Mazaheri@NASA.gov. 
t  Senior  Research  Scientist,  110  Exploration  Way,  Hampton,  VA  23666. 
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The  resulting  scheme  was  applied  to  various  time-dependent  problems,  demonstrating  second-order  accuracy 
for  the  solution  and  the  gradient  achieved  at  all  interior  and  boundary  nodes  in  uniform  and  nonuniform 
grids  at  every  physical  time  step,  and  rapid  convergence  for  solving  implicit-residual  equations  by  Newton’s 
method  (i.e.,  less  than  5  iterations  per  physical  time  step),  which  is  possible  by  the  compactness  of  the  RD 
schemes.  As  a  consequence  of  the  first-order  re-formulation  of  the  equation,  the  number  of  linear  relaxations 
performed  at  every  Newton  iteration  was  shown  to  increase  only  linearly  with  the  grid  size,  not  quadratically 
as  typical  for  diffusion  problems.  The  efficiency  of  the  developed  second-order  schemes  was  demonstrated 
for  linear  and  nonlinear  advection-diffusion  problems  on  highly  refined  grids,  up  to  30000  nodes. 

In  this  paper,  we  propose  a  very  simple  extension  of  the  second-order  schemes  to  higher-order.  We 
show  that  high-order  spatial  accuracy  can  be  achieved  simply  by  modifying  the  source  term  discretization. 
There  are  two  approaches  to  the  source  term  discretization:  1)  reformulation  of  the  source  terms  with  their 
divergence  forms,  and  2)  correction  to  the  trapezoidal  rule  for  the  source  term  discretization.  The  former 
technique  is  based  on  the  divergence  formulation  of  source  terms  proposed  in  Ref.:  write  the  source  term 
in  the  divergence  form  and  discretize  it  in  the  same  way  as  the  flux  divergence  term.  The  latter  is  based 
on  a  high-order  correction  to  the  trapezoidal  rule,  and  thus  called  here  the  generalized  trapezoidal  rule. 
In  either  case,  high-order  accuracy  is  achieved  by  making  low-order  truncation  error  terms  proportional  to 
the  residual,  which  thus  vanish  in  the  steady  state  and  yield  high-order  accuracy.  We  solve  the  resulting 
implicit-residual  equations  by  an  implicit  solver  based  on  the  second-order  Jacobian  matrix  developed  in  the 
previous  work.  As  we  will  show,  the  implicit  solver  is  as  powerful  as  Newton’s  method;  e.g.,  eight  orders 
of  magnitude  reduction  can  be  achieved  in  10  iterations.  To  enable  time-accurate  computations,  we  employ 
high-order  versions  of  the  backward  difference  formulas  (BDF),  which  are  incorporated  as  source  terms,  and 
solve  the  implicit-residual  equations  by  the  implicit  solver  over  each  physical  time  step.  In  this  manner,  the 
steady  state  is  made  equivalent  to  the  next  physical  time  with  all  the  benefits  of  the  hyperbolic  method 
retained. 

The  high-order  RD  schemes  developed  in  this  work  are  significantly  different  from  other  high-order  RD 
schemes  in  that  our  schemes  are  based  on  the  first-order  hyperbolic  system  formulation  of  the  advection- 
diffusion  equation.  In  this  approach,  the  loss  of  high-order  accuracy  in  the  intermediate  Reynolds  number, 
as  discussed  in  Refs.,'1-  cannot  occur  because  the  advective  and  diffusive  terms  are  fully  integrated  into  a 
single  hyperbolic  system.  If  the  original  advection-diffusion  equation  is  discretized,  a  high-order  RD  scheme 
needs  to  be  developed  for  the  diffusion  term  (i.e.,  second  derivative)  and  then  carefully  combined  with  an 
advection  scheme,  e.g.  by  using  a  blending  parameter  as  described  in  Ref.,  to  avoid  the  loss  of  accuracy. 
Furthermore,  while  high-order  RD  schemes  based  on  high-order  elements  require  extra  degrees  of  freedom  for 
each  variable,  our  schemes  are  based  on  linear  elements  for  any  order  of  accuracy  but  require  extra  gradient 
variable  to  be  added  to  the  solution  vector.  Note  that  the  number  of  extra  variables  in  the  high-order  elements 
increase  for  higher-order  accuracy,  but  the  number  of  extra  variables  required  in  our  approach  is  fixed  and 
independent  of  the  order  of  accuracy.  Our  approach  is  somewhat  similar  to  those  in  Refs.,1J-  but  again 
is  significantly  different  by  the  use  of  first-order  hyperbolic  system  formulation  of  the  advection-diffusion 
equation  and  by  the  source  term  discretization  techniques.  It  is  emphasized  that  our  schemes  require  only 
the  first  and  second  derivatives  of  the  source  term,  or  in  some  cases  the  first  derivatives  only;  they  do  not 
require  the  gradient  computation  for  the  solution  variables. 

The  third-order  schemes  developed  in  this  paper  are  similar  to  the  third-order  finite-volume  scheme  of 
Katz  and  Sankaran  ■  in  that  the  second-order  truncation  error  is  eliminated  by  making  it  proportional  to 
the  residual  and  the  upgrade  is  achieved  by  second-order  accurate  gradients.  However,  as  we  demonstrate 
in  this  paper,  the  proposed  high-order  RD  schemes  have  several  superior  features:  1)  implicit  solver  can  be 
constructed  by  the  Jacobian  derived  from  a  compact  second-order  RD  scheme,  2)  gradient  computations 
are  required  for  the  source  terms  only,  and  not  for  the  solution,  3)  stiffness  due  to  the  second  derivative  of 
the  diffusion  term  is  completely  eliminated,  4)  higher-order  schemes  can  be  constructed  beyond  third-order, 
5)  the  same  order  of  accuracy  is  achieved  for  the  gradients,  as  well.  In  particular,  the  fourth-order  scheme 
constructed  in  Section  5  is  remarkably  more  efficient  because  it  does  not  require  second  derivatives  of  the 
source  term,  which  are  required  in  the  schemes  described  in  Refs.  ■ 

In  this  paper,  we  focus  on  one-dimensional  linear  and  nonlinear  advection-diffusion  problems.  It  cer¬ 
tainly  serves  as  a  basis  for  the  development  of  high-order  multi-dimensional  RD  schemes  for  more  complex 
equations.  Yet,  more  importantly,  the  one-dimensional  high-order  schemes  developed  in  this  paper  could  po¬ 
tentially  bring  significant  improvements  to  practical  problems  such  as  material  thermal  response  calculations 
of  thermal  protection  systems  of  atmospheric  entry  vehicles,1'-1  and  the  experimental  aeroheating  data  re- 
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duction,-  ’  1  which  are  based  on  one-dimensional  analyses  and  routinely  used  in  industries  (e.g.  NASA).  The 
extension  to  higher  dimensions  is  beyond  the  scope  of  the  paper;  it  will  be  addressed  in  a  subsequent  paper. 

The  paper  is  organized  as  follows.  In  the  next  section,  the  time-dependent  hyperbolic  advection-diffusion 
system  is  described.  In  Section  3,  a  compact  second-order  residual-distribution  scheme,  a  steady  solver,  and 
the  second-order  discretization  are  discussed.  In  Section  4,  the  third-  and  fourth-order  RD  schemes  with 
source  term  reformulation  are  proposed.  In  Section  5,  the  third-,  fourth,  and  sixth-order  RD  schemes  with 
source  term  discretization  are  developed  and  proposed.  Numerical  results  are  then  presented  in  Section  6. 
Finally,  Section  7  concludes  the  study  with  remarks. 


II.  Time-dependent  Hyperbolic  Advection-Diffusion  System 


We  start  with  a  linear  advection-diffusion  equation  to  simplify  the  discussion.  We  will  extend  the  discus¬ 
sion  later  to  a  more  general  nonlinear  advection-diffusion  equation. 

Consider  the  one-dimensional  (1-D)  time-dependent  advection-diffusion  equation: 

dtu  +  adxu  =  isdxxu  +  S(x),  (1) 

where  a  and  is  are  both  taken  to  be  positive  constant,  and  S  is  the  source  term.  We  will  follow  the  procedure 
we  described  in  Ref.  and  rewrite  the  above  equation  as  a  first-order  hyperbolic  advection-diffusion  system: 

a 

dTu  =  —  a  dxu  +  is  dxp  —  —  u  +  S(x),  (2) 

dTp  =  ( dxu-p)/Tr ,  (3) 


where  the  relaxation  time,  Tr  >  0,  is  arbitrary  and  defined  as  described  in  Ref.,  and  S  includes  any  existing 
source  terms  from  the  advection-diffusion  problem,  S,  as  well  as  any  additional  terms  that  arise  from  the 
implicit  time-stepping  scheme,  At  is  the  physical  time  steps,  and  r  is  the  pseudo  time  step.  Note  that  the 
dtp  is  taken  as  pseudo  time  derivative,  dTp. 

The  variable  a  depends  on  the  order  of  the  Backward-Differencing-Formula  (BDF);  1  for  the  lst-order 
(BDF1),  3/2  for  the  second-order  (BDF2),  11/6  for  the  third-order  (BDF3),  25/12  for  the  fourth-order, 
and  147/60  for  the  sixth-order  time  discretizations.  The  remaining  terms  in  the  BDF  are  stored  in  the 
source  term  function  S{x).  Towards  the  pseudo  steady  state  state,  the  variable  p  approaches  the  solution 
gradient  and  hence  the  above  equation  becomes  identical  to  the  original  advection-diffusion  equation  with 
the  time  derivative  discretized  by  the  BDF,  i.e.,  a  consistent  discretization  of  the  original  equation.  The 
system  reduces,  of  course,  to  the  original  steady  equation  also  in  the  physical  steady  state  when  dtu  =  0. 
Second-order  accurate  unsteady  computations  have  been  demonstrated  based  on  the  above  formulation  in 
Ref.6 

In  the  vector  form,  our  time-dependent  first-order  advection-diffusion  system  can  be  written  as 


where 


<9U 

Ifr 


<9U 

L  dx 


=  S, 


(4) 


u  = 

U 

,  A  = 

a  —v 

,  S  = 

— au/At  +  S(x) 

V 

—  1/Tr  0 

-p/Tr 

For  any  positive  Tr,  the  Jacobian  has  the  following  two  real  eigenvales: 


Ai  = 


1-W1  + 


4v 

a?Tr 


Aa  — 


1  +  1/1 


4  v 
a^fr 


(5) 

(6) 


with  linearly  independent  eigenvectors,  and  therefore  the  above  system  is  hyperbolic  in  the  pseudo  time. 
The  hyperbolicity  serves  mainly  as  a  guide  for  discretization:  various  discretization  techniques  are  available 
for  hyperbolic  systems,  e.g.,  upwinding.  In  addition  to  the  convenience  in  discretization,  the  major  benefits 
are:  1)  the  hyperbolic  discretization  results  in  a  strong  coupling  between  the  two  variables  that  results  in  the 
equal  order  of  accuracy  for  both  u  and  p  =  dxu  on  arbitrary  grids,  and  2)  0(l/h)  acceleration  is  achieved 
in  iterative  convergence  for  the  linearized  residual  equation  in  implicit  solvers  due  to  the  0(l/h)  condition 
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number  of  the  Jacobian,  which  is  0(h)  smaller  than  that  of  the  Jacobian  derived  from  a  conventional 
diffusion  scheme.  For  explicit  schemes,  the  0(1 /h)  acceleration  is  achieved  by  the  0(h)  time  step  typical  in 
methods  for  hyperbolic  systems,  which  is  0(1 /h)  times  larger  than  a  typical  time  step  of  0(h 2)  for  diffusion. 
The  0(l/h)  acceleration  in  convergence  has  been  demonstrated  over  traditional  methods  for  the  diffusion 
equation, !  ’  for  the  advection-diffusion  equation, 2’  for  the  compressible  Navier-Stokes  equations,  and  most 
recently  for  time-dependent  linear  and  nonlinear  advection-diffusion  problems  using  RD  technique. 

In  the  next  sections,  we  first  briefly  describe  the  second-order  time-dependent  discretization  process  using 
the  RD  scheme.  Further  details  are  given  in  Ref.1'  We  then  extend  the  order  of  accuracy  of  the  scheme  to 
third-order  and  fourth-order  RD  hyperbolic  advection-diffusion  schemes  with  reformulation  of  the  hyperbolic 
system,  and  finally  to  higher  order  (fourth-  and  sixth-order)  RD  hyperbolic  advection-diffusion  schemes  with 
introduction  of  new  source  term  discretization  technique. 


III.  Second-Order  RD  Hyperbolic  Advection-Diffusion 

The  RD  method  requires  1)  evaluation  of  the  cell  (or  element)  residuals,  and  2)  the  distribution  of  the 
residuals  to  the  nodes  bounding  the  cell.  Consider  a  one-dimensional  domain  discretized  with  N  nodes 
that  are  distributed  arbitrary  over  the  domain  of  interest  with  the  solution,  u,  and  the  solution  gradient, 
p  =  dxii,  data  stored  at  the  nodes  denoted  by  a *=1,2, 3, ...,  N  (Fig.l).  We  define  the  cell-residual  <&c  by 

hL  hR 

- o - o - o - 

i-1  i  i+1 


Figure  1.  Schematic  of  grid  spacing  for  a  1-D  grid. 


integrating  the  spatial  part  of  the  Eq.  (4)  over  the  cell,  C ,  defined  by  the  nodes  *  and  *  +  1: 


(-AU*  +  S  )dx, 


(7) 


{~a(ui+ 1  -  Ui)  +  u(pi+ 1  -  p^ 
Y  [ui+ 1  ~  «*)] 


(8) 


The  first  term  of  the  Eq.  (8)  is  the  result  of  the  exact  integration  of  -AUj,  over  the  cell  C.  The  integration 
of  the  second  term  is  not  exact  and  therefore  is  the  source  of  the  overall  discretization  error.  We  will  further 
discuss  this  discretization  error  in  the  following  sections. 

We  construct  the  spatial  discretization  of  Eq.  (4)  by  distributing  the  cell  residuals  to  the  nodes: 


^  = -f(S+$L +  £"$*),  (9) 

clt  hi 

where  <&L  and  denote  the  cell-residuals  over  the  left  and  right  cells  of  the  node  i,  respectively,  and  hi  is 
the  dual  volume  (see  Fig.  1)  defined  by 


hi  = 


hL  +  hR 
2 


hL  %i  %i—  I?  hji  —  1 


(10) 


Note  that  the  pseudo  time  derivative  on  the  left  hand  side  is  retained  here  just  for  the  sake  of  illustration.  It 
will  be  ignored  in  the  implicit  formulation  that  follows.  Our  aim  is  to  directly  solve  it  for  the  pseudo  steady 
state,  which  corresponds  to  the  next  physical  time.  In  this  sense,  our  method  is  not  a  dual-time  stepping 
method.  The  distribution  matrices  B~  and  B+  (Fig.  2)  typically  do  not  affect  the  order  of  accuracy  of 
the  discretization  scheme;  it  is  determined  by  the  cell-residuals, and  therefore  our  focus  will  be  on  the 
residual  evaluation.  (Refer  to  Ref.1  for  formulation  of  the  distribution  matrices  B~  and  B+). 
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Figure  2.  Residual  distribution  to  the  nodes. 


A.  Second-Order  Discretization 

We  may  use  a  simple  trapezoidal  rule  for  the  integration  of  the  source  term  S  over  the  cell  C  and  arrive  at 
the  following  cell  residuals: 


<j>c 


-a(ui+ 1  -  Ui)  +  v(pi+ 1  -  Pi)  -  (xi+i 


•Ei) 


a 
A t 


(tii+l  +  'U'i)  /  2 


fe+1 


(^2+1  Ui  (Xi- |_i  (.Pi+l  +  Pi) /2) 


"I-  ^i)/^ 

0 


n—  l,n 


(ii) 


where  k  and  n  are  the  Newton  iteration  counter  (as  described  below)  and  the  physical  time  index,  respectively. 
Note  that  the  second  term  of  the  Eq.  (11),  which  is  computed  at  the  two  previous  physical  time  steps,  is 
constant  during  the  Newton  iteration,  and  thus  will  not  contribute  to  the  Jacobian. 

The  implicit  formulation  is  defined  by 


U fc+1  =  Ufc  +  AUfc, 


(12) 


where  U  =  (ui,pi,v,2,P2,  ■  ■  ■  ,un,Pn )  and  k  is  the  iteration  counter.  The  correction  AUfc  =  Ufc+1  —  Ufc  is 
determined  as  the  solution  to  the  linear  system: 


5ResAUfe 

dU 


— Res'b 


(13) 


where  Resfc  is  the  right  hand  side  of  Equation  (9),  which  is  the  unsteady  residual  vector  evaluated  by  Ufc. 
Note  that  the  pseudo  time  derivative  has  been  ignored.  The  Jacobian  matrix  is  sparse,  having  three  2x2 
blocks  in  each  row  for  all  interior  nodes  and  two  blocks  for  boundary  nodes.  The  i-th  pair  of  rows  of  the 
linear  system  is  given  by 

Ji-iAU-L1  +  JiAU f  +  Ji+iAUf+1  =  -1(S+$L  +  £T$*)fe,  (14) 

fli 


where  A  U?_!  =  (Au^AplJ,  A  U?  =  (Au^  Atf),  A  U?+1  =  (Au,*+1,  Ap*+1), 

1  d(B+<5>L ) 
i_1  “  hi  dUi-i  ’ 

1  (d(B+<S>L)  {dB~<f>R)\ 

1  ~  hi{  dUi  +  dUi  )  ’ 

1  d(B~<i>R) 

1+1  hi  dUi+1 


(15) 

(16) 
(17) 


We  note  that  the  derivative  of  the  distribution  matrices,  B~  and  1?+,  are  zero  for  linear  advection-diffusion 
problems  and  for  nonlinear  problems  with  constant  a/v  values,  which,  for  simplicity,  are  the  only  nonlin¬ 
ear  system  considered  in  this  paper.  However,  the  proposed  schemes  are  applicable  to  general  nonlinear 
advection-diffusion  problems. 
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The  residual  Jacobians  needed  in  the  Newton  solver  are  exactly  in  the  same  form  as  the  above  equations, 
but  the  derivatives  of  the  cell-residuals  now  include  the  contribution  from  the  physical  time  derivative: 


d®R 


d$L 

dUt 


a  —  Iir—— 

2  At 

—v 

1 

hR 

_  TV 

_2T ; 

,  a 

—a  —  nL— — 

2 At 

V 

1 

hL 

TV 

2  Tr 

(18) 


(19) 


At  each  physical  time  level  n,  we  solve  the  pseudo  steady  problem  by  Newton’s  method  with  the  current 
solution  at  n  as  the  initial  solution.  This  results  in  second-order  discretization,  which  was  demonstrated 
with  several  examples  (linear  and  nonlinear)  in  Ref/ 

We  now  focus  on  the  truncation  error  of  the  source  term  discretization  as  it  is  needed  for  our  discussion 
in  the  next  sections  where  we  show  two  different  approaches  to  obtain  higher-order  RD  schemes.  We 
demonstrate  the  truncation  error  of  the  system  by  considering  the  residual  of  the  second  equation  in  our 
hyperbolic  advection-diffusion  system: 


=  (ui+ 1  -  Ui  -  h{jpi+ 1  +pi)/2)/Tr, 


(20) 


where  h  is  the  width  of  the  cell  C  (i.e.  Xj+i  —  xi).  We  now  expand  the  around  the  node  i  to  obtain  the 
truncation  error  ( T.E .)  for  the  second  equation,  9rp,  which  after  some  algebra  and  rearrangements  becomes: 


T.E.  (dTp)  = 

=  o(h2) 


h  h 2 

( dxUi  Pi)  T  {dxxUi  dxPi)  T  (A 

2  o 


xxx'Ui  ^&xxPi)  d-  0(h  ) 


/Tr 


(21) 


where  the  first  two  terms  will  simply  vanish  as  they  satisfy  the  original  equation  in  the  pseudo  steady  state. 
The  non-vanishing  last  term  makes  the  current  discretization  second-order.  The  same  results  are  obtained 
with  the  first  equation  and  therefore  is  not  repeated  here.  The  above  analysis  shows  that  the  error  arises 
from  the  source  term  discretization,  not  from  the  flux  balance  term,  which  is  exact  in  one  dimension.  It 
implies  that  we  can  achieve  high-order  accuracy  only  by  improving  the  source  term  discretization.  In  fact, 
in  Ref.,1  :  a  fourth-order  finite-difference  scheme  is  constructed  in  this  manner,  i.e.,  by  a  high-order  source 
term  discretization  with  explicit  pseudo-time  stepping  targeted  for  steady  problems. 

Note  that  the  above  truncation  error  analysis  is  based  on  the  cell-residual  expanded  at  a  node.  Hence, 
the  nodal  residual  has  the  same  order  as  the  cell-residuals  because  it  is  constructed  as  a  weighted  sum  of  the 
cell- residuals  as  shown  in  Eq.  (9).  Note  also  that  the  truncation  error  analysis  based  on  the  cell-residual  is 
local  to  the  cell  and  thus  valid  for  any  size  of  the  cell,  i.e.,  valid  for  uniform  and  nonuniform  grids.  In  the 
rest  of  the  paper,  the  truncation  error  analysis  is  all  based  on  the  cell-residual. 

In  the  next  sections,  we  show  two  different  techniques  that  will  lead  to  higher-order  schemes.  In  partic¬ 
ular,  we  discuss  on  1)  reformulation  of  the  source  term  with  its  divergence  form,  and  2)  new  source  term 
discretization  technique  which  acts  as  a  correction  to  the  trapezoidal  rule. 


IV.  High-Order  RD  Scheme  with  Source  Term  Reformulation  (Third-  &; 

Fourth-Order) 

Consider  the  source  term,  S ,  as  the  divergence  of  a  function  fs:  =  S.  We  now  replace  the  source 

terms  with  their  divergence  forms,  and  rewrite  our  first-order  hyperbolic  advection-diffusion  equation  as 

dTu  =  -adxu  +  vdxp  +  dxfS,  (22) 

dTp  =  (dxu  —  dxfp  )/Tr,  (23) 

where  and  are  the  divergence  formulation  of  the  source  terms  in  the  u  and  the  p  equations,  respectively. 
With  the  above  reformulation  of  the  advection-diffusion  equation,  which  will  be  discussed  in  more  details, 
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the  residual  evaluations  of  the  system  becomes  exact  with  no  special  discretization  scheme: 


(-AVx  +  fsx)dx, 


(24) 


!-a{ui+ 1  -  Ui)  +  v(pi+i  -  Pi)  +  (fu)i+ 1  -  {fu)i 

1  •  (25) 

—  [ui+i  - Ui  -  (fp)i+ 1  +  {fp)i)] 

We  remark  that  even  though  the  residual  evaluation  of  our  reformulated  advection-diffusion  system  is  exact, 
the  overall  accuracy  of  the  scheme  depends  on  the  accuracy  of  the  divergence  formulation  of  the  source  terms 
and  how  it  is  formulated. 


A.  Third-Order  Scheme  with  Divergence  Form 

Following  the  divergence  formulation  introduced  in,  we  write  the  source  flux  in  a  more  general  form: 


fS 


m— 3  ( _i  'in— 1 

V  1  \  (x  —  Xi)ndxn-lS, 

Z '  71. 1 


n—1 


(26) 


=  (x  ~  Xi)S -  Xi)2dxS +^(x  -  Xi)3dxxS  -  -  Xi)4dxxxS  +  . ..,  (27) 

where  the  source  term  S  and  its  derivatives  Sx,  Sxx,  etc.  are  not  constants  but  functions  that  vary  in 
space.  We  remark  that  the  source  flux,  fs,  is  not  necessarily  a  polynomial  of  order  m  because  it  depends 
on  the  derivatives  of  the  source  term  S.  In  this  paper,  we  do  not  consider  high-order  schemes  that  require 
evaluation  of  third  or  higher  derivatives.  The  divergence  of  the  source  flux,  dxfs ,  for  the  third-order  accuracy 
(i.e.  m  =  3)  is 

dxfs  =  S+^-dxxxS  =  S  +  0(h3),  (28) 

o 


which  is  identical  to  the  original  equation  up  to  the  third-order.  We  remark  that  this  error  does  not 
necessarily  limit  the  maximum  order  attained  by  numerical  schemes.  In  the  next  section,  we  will  show  that 
a  fourth-order  scheme  can  be  constructed  by  a  slightly  modified  divergence  form  with  m  =  3. 

Discretization  of  the  reformulated  hyperbolic  system  leads  to  the  cell  residual,  for  example,  for  the  second 
equation  as 


=  [ui+i-Ui-{fp)i+i  +  (fp)i]/Tr, 


(29) 


where  special  care  must  be  taken  when  we  evaluate  ff  and  ff,:i  (because  of  the  presence  of  aq  in  the 
divergence  formulation)  as  they  depend  on  whether  the  cell  residual  is  being  distributed  to  the  left  or  to  the 
right  node  of  the  cell  (Fig.  3): 


[«i+l  “  Ui  -  Up)i+ 1  +  (fp)i\  /Tr, 


h2  h3 

^2+1  aq  hftSi- )_i  +  3X  ~^~^xx^i+i 


/Tr, 


[Ui+ 1  -Ui  -  (fp)z+ 1  +  (fp  )i]  /Tr, 
h?  It ^ 

'Ui  ~~p~dxSi  +  ~z~^xx^i 

2  D 


/Tr. 


(30) 


(31) 


The  source  term  discretization  is,  therefore,  not  conservative,  which  is  natural  and  appropriate  because  the 
source  term  does  not  have  a  conservative  property  and  should  not  be  discretized  in  such  a  way.  If  it  is 
conservative,  the  global  sum  of  the  cell-residual  for  the  source  term  will  depend  only  on  the  boundary  data 
(telescoping  property),  which  is  wrong  for  the  source  term. 

We  now  show  that  the  truncation  error  ( T.E .)  of  the  resulting  RD  scheme  with  the  above  divergence 
formulation  of  the  source  term  (with  to  =  3)  is  in  fact  third-order.  Again,  for  demonstration  purposes  we 
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Figure  3.  Cell  residual  evaluation  with  divergence  form  of  the  source  term. 


consider  the  second  equation  (i.e.  dTp)\  the  same  process  can  be  repeated  for  the  first  equation.  We  first 
expand  the  source  fluxes  around  the  node  i: 


j  dxfs dx  =  ff+  r  -  f?  =  hRdxfs  +  +  ^d^f3  +  0(h4) 


where  the  dxfs  is  given  by  Eq.  (28),  and 


d~'s  =  «.s  +  fa*.s+fa—s+o(/.‘), 

S  2  h3  4 

8xxxf  =  8XXS  -t-  hRdxxxS  T  hRdxxxxS  T  —^~dxxxxxS  T  0{h  ). 


(32) 

(33) 

(34) 


Using  Eq.  (29),  we  evaluate  the  truncation  error  of  the  equation  dTp  by  substituting  Eqs.  (33)  and  (34)  into 
Eq.  (32)  and  expanding  all  the  terms  in  Eq.  (29)  around  the  node  i: 


T.E.  (dTp)  = 


{dxUi  Pi)  n  {^xxU'i  dxPi)  +  (dxxxUi  dxxPi) 
2  o 


/Tr 


h3 

— rr(8xxxxUi  14 dxxxPi)  "b  0(h  ) 


24 

=  0  +  0(h3) 


/Tr 


(35) 


where  the  first  three  terms  vanish  as  they  satisfy  our  original  equation.  The  presence  of  the  last  term  confirms 
that  the  proposed  divergence  formula  with  m  =  3  makes  our  scheme  third-order  accurate.  The  same  result 
is  obtained  for  the  first  equation  with  the  divergence  formulation  of  the  corresponding  source  terms,  and 
therefore  the  process  is  not  repeated  here. 

We  remark  that  the  cost  of  this  new  third-order  accurate  RD  hyperbolic  advection-diffusion  scheme  using 
the  divergence  formulation  of  the  source  terms  is  only  the  evaluation  of  the  first  and  second  derivatives  of 
the  source  terms;  i.e.  dxS  and  dxxS.  Note  that  the  source  flux  is  third-order  accurate  (see  Eq.  28).  Thus, 
according  to  the  Eqs.  (30)  and  (31),  we  need  second-order  and  first-order  accurate  discretization  for  dxS  and 
dxxS ,  respectively.  We  derive  these  discretization  by  constructing  a  quadratic  function  between  the  nodes 
i  and  its  two  neighbors  i  —  1  and  i  +  1  to  arrive  (after  some  algebra)  with  the  following  formulas  that  are 
applicable  to  the  internal  nodes  of  general  arbitrary  grids  (uniform  and  nonuniform): 


dxSi  =  '12rSi~1  ,+  kL]S\+  kLSl+1  +  Q(h2) 


(d XX  — 


hRhL{hR  +  hL) 

hRSj-i  —  (hR  +  hi,)Si  +  h.LSj+1 
hRh.L(hR  +  hL)/2 


■0(h), 


(36) 

(37) 


where  hR  and  hR  are  defined  in  Eq.  (10).  The  corresponding  one-sided  formulas  for  the  boundary  nodes  of 
general  arbitrary  grids  are 


dxS1 

8xSn 

8Xx^i 

8xxSn 


-hL.(hL.  +  2 hR)S1  +  (hR  +  hL*)2S2  -  h2RS3 


hRhL*(hR  +  hL .) 
hL(h-L  +  2hR*)SN  —  (hR*  +  Jir)2Sn-i  +  hR,S]y- 1 
hR-hL(hR •  +  hL) 

h,L*S\  —  (hR  +  hR*)S2  +  hRS3  n(u\ 
hRhL.(hR  +  hL.)/2  +  1 

—hR*S]y  +  (hR*  +  /i_l)5jv-i  —  2 

hR*h,L(hR *  +  hi)/ 2 


0(h2), 


o(h2), 


0(h), 


(38) 

(39) 

(40) 

(41) 
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where  hi*  and  Hr*  are  defined  as 


hi,*  =x3-  x2,  hR*  =  xN-i  -  xN-2-  (42) 

It  is  clear  from  the  above  formulas  that  each  derivative  can  be  computed  in  a  three-point  stencil.  Con¬ 
sequently,  the  residual  at  a  node  is  defined  in  a  five-point  stencil  in  the  interior,  a  four-point  stencil  at  the 
nodes  adjacent  to  the  boundary,  and  a  three-point  stencil  at  the  boundary  nodes.  In  the  next  section,  we 
demonstrate  that  a  fourth-order  scheme  can  be  constructed  without  extending  the  stencil. 


B.  Fourth-Order  Scheme  with  Divergence  Form 

Here,  we  show  how  a  simple  modification  to  the  presented  divergence  form  upgrades  the  scheme  order  by 
one  order;  that  is  our  third-order  scheme  becomes  fourth-order  with  no  additional  cost.  To  gain  an  order, 
we  propose  the  divergence  form  of  the  source  terms  (i.e.  function  fs ),  to  be  defined  as 


m>  3 


/s  =  E 


(-i) 


n—  1 


(x-  x)ndxn-iS, 

1.  „  1,  „  1 


=  (x-  x)S  -  ^(x  -  x)2dxS  A  -(x  -  x)3dxxS  -  —  (x  -  x)AdxxxS  A  . . . , 


(43) 

(44) 


where  x  =  (x^  A  Xj+i)/2.  Note  that  the  previous  divergence  form  was  formed  around  the  x,:  while  this 
new  formulation  defines  the  divergence  function  around  the  mid-point  x.  Similar  to  the  divergence  form 
presented  in  the  previous  section,  the  source  term  and  its  derivatives  are  not  constants  but  functions  that 
vary  in  space.  And  therefore,  when  the  source  flux  is  differentiated  it  recovers  the  original  source  term  up 
to  0(hm)  around  the  node  i.  For  m  =  3,  we  have: 


dxfs  =  S+{x  RX)3dxxxS  =  S  +  0(h3).  (45) 

o 

Again,  it  does  not  limit  the  maximum  order  of  accuracy  for  numerical  schemes,  and  it  is  possible  to  construct 
even  higher-order  schemes. 

We  now  show  that  this  third-order  source  divergence  formulation  will  result  in  a  fourth-order  accurate 
RD  scheme.  We  do  this  similarly  to  the  process  explained  in  the  previous  section  except  that  the  second 
and  third  derivatives  of  the  fs  are  now  defined  as: 


dxxfs  =  dxs+{x  2x)2 dxxxs+{x  6x)3dxxxxs  +  o(h4), 


&xxxf  —  &XXS  A  (x  x)dxxxS  A  (x  x)  dxxxxS  A 


(x  —  x)c 

6 


dxxxxxS  A  0(h4). 


(46) 

(47) 


Again,  for  discussion  purposes,  we  consider  the  second  equation  of  the  hyperbolic  advection-diffusion  system, 
dTp.  The  truncation  error  of  the  p  equation  after  substituting  Eqs.  (46)  and  (47)  into  Eq.  (32)  becomes: 


T.E.  (dTp)  = 


h  h 2 

{dxUi  Pi )  A  ~R~{dxxUi  9xPi )  A  ~^~{OxxxUi  dxxPi ) 
2  6 


/Tr 


h3  hA  5 

-^^(^xxxx'U’i  & xxxPi )  4”  i2Q  iPxxxxx^i  xxxxPi )  4”  0(h  )  /Tr 


=  0  A  0(h4), 


(48) 


which  is  fourth-order  accurate  because  of  the  cancellation  of  the  first  four  terms  of  the  truncation  error 
equation.  Another  great  property  of  this  new  divergence  formulation  is  the  equal  evaluation  of  the  ff+1  —  f f 
regardless  of  whether  the  source  flux  is  being  transferred  to  the  node  i  or  i  A  1.  This  greatly  simplifies  the 
implementation  of  this  form  of  divergence  formulation. 

Note  that  the  identical  residual  is  distributed  to  the  left  and  right  nodes  in  this  scheme.  It  then  appears 
that  the  source  term  discretization  is  conservative.  However,  it  is  actually  not  conservative  for  the  source 
term  because  the  cell-residual  for  the  source  term  depends  on  the  coordinate  of  the  mid-point  of  the  cell 
and  thus  it  is  not  telescoping.  Again,  this  is  natural  and  appropriate  for  the  source  term,  which  has  no 
conservation  property  and  should  be  local. 
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V.  Higher-Order  RD  Scheme  with  Generalized  Trapezoidal  Rule  (Third-, 

Fourth-  &c  Sixth-order) 


In  the  previous  section,  we  reformulated  the  original  hyperbolic  advection-diffusion  system  with  a  gener¬ 
alized  divergence  form  of  the  source  terms  and  arrived  at  third-  and  fourth-order  RD  schemes.  We  specifically 
showed  that  both  the  third-  and  the  fourth-order  RD  schemes  developed  with  the  divergence  formulation  of 
the  source  terms  require  the  evaluation  of  the  first  and  second  derivatives  of  the  source  terms.  Fifth-  and 
higher-order  schemes  can  be  systematically  constructed  with  the  divergence  formulation,  but  may  require 
the  evaluation  of  third-  and  higher-order  derivatives  that  would  extend  the  stencil  to  the  neighbors  of  the 
neighbors  and  beyond.  From  a  practical  point  of  view,  such  high-order  schemes  are  not  very  attractive,  and 
thus  not  considered  here. 

In  this  section,  we  introduce  a  different  technique  to  develop  even  higher-order  schemes  without  the 
need  to  evaluate  gradients  beyond  the  second-derivatives.  We  also  show  that  with  this  new  technique  the 
fourth-order  RD  scheme  is  even  less  computationally  intensive  than  the  third-  and  fourth-order  RD  schemes 
that  are  developed  with  the  divergence  formulation  of  the  source  terms. 

Consider  the  vector  form  of  our  first-order  hyperbolic  advection-diffusion  system: 


<9U 

Ih 


<9U 
L  dx 


=  S. 


(49) 


We  showed  in  Section  III  that  the  source  term  discretization  with  the  trapezoidal  rule  provides  a  second-order 
accurate  scheme.  This  trapezoidal  rule  can  be  written  as 


Sdx  =  -f(SL  +  Sr), 


(50) 


where  for  the  second-order  scheme  Sl  =  Si  and  Sr  =  ;  i.e.  the  arithmetic  averaging  of  the  source  terms 

between  the  left  and  the  right  nodes  (Fig.  4). 


Figure  4.  Source  term  discretization. 

Generalizing  the  trapezoidal  rule,  we  propose  the  following  formula  for  the  left  and  right  source  terms, 
Sl  and  Sr  respectively: 

SL  -  Si  +  Ct'dxSi  +  C£dxxSi,  (51) 

SR  =  Si+1+C?dxSi+1+C?dxxSi+1 ,  (52) 

where  C\J  and  are  constants  of  0(h),  and  Cf  and  C 2  are  of  0(h2).  A  somewhat  similar  approach  is 
taken  in  Ref.,  introduced  for  upgrading  the  order  of  accuracy  of  a  finite- volume  scheme,  where  not  only  the 
source  term  but  also  the  interface  flux  need  to  be  modified  for  high-order  accuracy.  We  find  these  constants 
by  making  sure  that  the  nodal  residuals  of  the  two  equations  in  our  hyperbolic  system  are  accurate  up  to  the 
desired  order  of  accuracy  of  the  scheme  we  are  seeking  to  develop.  With  the  above  equation,  the  maximum 
possible  order  of  accuracy  is  six.  We  also  note  that  the  fifth-order  RD  scheme  with  the  above  proposed 
source  term  discretization  becomes  mathematically  sixth-order. 

A.  Third-Order  RD  Scheme  with  Generalized  Trapezoidal  Rule 

A  third-order  RD  scheme  can  be  obtained  if  the  coefficients  of  the  proposed  discretization  formulation  satisfy 
the  following  equations: 

Cf  +  Cf  =  0,  C^hR  +  C^  +  C^  =  -^,  C?hR  +  (53) 

which  can  be  derived  from  the  truncation  error  analysis.  Here,  there  are  four  unknowns  for  three  equations; 
thus  there  are  infinite  combinations  with  the  C ^  as  a  free  parameter.  We  also  note  that  there  is  no  relation 
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between  this  third-order  RD  scheme  and  the  third-order  scheme  developed  with  the  divergence  formulation; 
the  nodal  residual  of  the  source  terms  obtained  with  the  previous  third-order  scheme  depends  only  on  the 
information  from  its  immediate  node  neighbor  and  not  the  node  itself.  On  the  other  hand,  the  new  third- 
order  RD  scheme  that  is  developed  with  a  correction  to  the  trapezoidal  rule  depends  both  on  its  immediate 
node  neighbor  and  the  node  itself.  Consequently,  it  is  not  possible  to  reproduce  the  previous  third-order 
scheme  by  any  choice  of  the  coefficients.  Note,  however,  that  this  scheme  has  the  same  stencil  with  the 
previous  third-order  scheme:  a  five-point  stencil  in  the  interior,  a  four-point  stencil  at  the  nodes  adjacent 
to  the  boundary,  and  a  three-point  stencil  at  the  boundary  nodes.  The  same  is  true  for  the  fourth-order 
scheme  derived  in  the  next  section. 

Assuming  C ^  =  —Hr/ 6,  the  left  and  right  source  functions  of  the  source  discretization  may  be  chosen 
as: 

SL  =  Si+bfdxSi^+^d^Si  (54) 

SR  =  Sl+1  -  ^dxSi+1  -  ^-dxxSl+1.  (55) 

o  1000 

We  can  now  evaluate  the  truncation  error  of  the  hyperbolic  advection-diffusion  system.  Here  we  show  the 
truncation  error  of  the  second  equation;  the  same  will  be  true  for  the  first  equation  as  well.  The  cell  residual 
of  the  second  equation  with  the  third-order  RD  scheme  with  the  generalized  trapezoidal  rule  is: 

=  ~a(ui+1  -  Ui)  +  u(pi+1  -  pi)  +  ~y(Sl  +  Sr), 

=  -a(ui+ 1  -  ui)  +  v{pi+ 1  -  Pi)  +  +  Si) 

-  {dxSi+1-dxSi)--^{dxxSi+1-dxxSi).  (56) 

We  then  expand  the  cell  residual  around  the  node  i  to  arrive  with  the  following  truncation  error: 

T.E.(drUi)  —  (  adxUi  H-  vdxpi  H-  S H-  ^  (  ^^xx^i  vdXxPi  &xSi) 

h2 

H-  — (  QJ&xxx'UJi  “I”  V^xxxPi  “1“  ^xx^i) 

6 

h 3 

+  xxxxPi  H-  0.988  &xxx  )  +  0(h 4) 

=  0  +  0(h3),  (57) 

where  the  first  three  terms  vanish  but  the  last  term  makes  the  scheme  third-order  accurate.  As  we  will  show 
next,  smaller  values  for  the  coefficients  C-f  or  move  the  schemes  to  fourth-order  accurate.  This  is  further 
shown  in  the  following  section. 

B.  Fourth- Order  RD  Scheme  with  Generalized  Trapezoidal  Rule 

Fourth-order  accuracy  is  achievable  by  a  simple  adjustment  to  the  condition  given  in  Eq.  (53): 

C?  +  C?  =  0,  C*hR  +  CL  +  C*  =  -hf ,  °thR  +  C*  =  -^ |.  (58) 

6  2  12 

Again,  the  coefficients  cannot  be  determined  uniquely,  and  there  are  many  solutions.  A  particularly  inter¬ 
esting  solution  is  the  following: 

cf  =  +^,  cf  =  0,  Cf  =  -^,  Cf  =  0.  (59) 

As  will  be  shown  shortly,  it  leads  to  a  fourth-order  RD  scheme  with  the  source  term  discretization  by  Eqs.  (51) 
and  (52): 

SL 

Sr 
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It  is  remarkable  that  fourth-order  accuracy  is  achieved  without  the  second  derivatives,  which  is  less  expensive 
than  the  third-order  schemes  developed  in  the  previous  sections.  It  is,  of  course,  possible  to  develop  a  fourth- 
order  RD  scheme  with  addition  of  the  second  derivatives  in  the  source  term  discretization.  For  example  the 
following  constants  will  also  result  in  a  fourth-order  RD  scheme: 


C 


L 

1 


ci 


24 


C 


R 


hR  r<R  _  ,  hR 
”T’  °2  ~  +  24’ 


(62) 


which  makes  this  fourth-order  RD  scheme  identical  to  the  fourth-order  RD  scheme  constructed  by  the 
divergence  formulation  in  Section  B.  Although  it  is  possible  to  construct,  these  fourth-order  schemes  are  not 
very  attractive  as  they  require  second  derivatives  of  the  source  term,  and  therefore  not  considered  in  this 
paper. 

We  now  prove  the  fourth-order  accuracy  of  the  scheme  by  evaluating  the  cell  residual  of,  for  example, 
the  first  equation;  i.e.  : 


=  ~a(ui+i  -  Ui)  +  u(pi+1  -  pi)  +  ~^{SL  +  Sr), 

h  h? 

=  -a(ui+1  -  Ui)  +  v(pi+1  -p^  +  -^(Si+1  +  Si)  -  ~^(dxSi+i  -  dxSi).  (63) 

Expanding  the  cell  residual  around  the  node  i.  we  obtain  the  truncation  error  of  the  first  equation  (after 
some  algebra)  as 


T .E .(dTUi)  —  (  ddxUi  -\~  lsdxPi  “t~  *5))  -t-  ^  (  d^xx^i  T  R^xxPi  d-  ^x^i) 
h2 

d"  ~ TT  (  0‘dxxx'U'i  d”  R^xxxPi  d”  &xx^i) 

6 

h3 

+  Q’dxxxx'Ri  d”  vdxxxxPi  d”  ^xxx^i) 

/l4  5 

d“  (  R^xxxxx'Ri  d"  R^xxxxxPi  d*  R^xxxx^i')  d”  0{}P  ) 


120 
0  +  O(h4), 


6 


(64) 


where  the  first  four  terms  of  the  above  equation  vanish  (because  of  consistency  relations).  Similar  conclusion 
is  obtained  for  the  second  equation  and  therefore  it  is  not  repeated  here.  Also,  note  that  the  accuracy 
is  achieved  for  general  arbitrary  grids.  The  additional  cost  of  upgrading  the  second-order  scheme  to  the 
fourth-order  RD  scheme  is  only  due  to  the  evaluation  of  the  second-order  accurate  first  derivative  of  the 
source  terms.  The  general  second-order  accurate  first  derivative  formulation  for  arbitrary  grids  is  provided 
in  Sec.  A. 


C.  Sixth-Order  RD  Scheme  with  Generalized  Trapezoidal  Rule 

In  this  section,  we  present  a  new  sixth-order  RD  scheme  with  relatively  similar  cost  to  our  newly  introduced 
third-order  RD  schemes.  Specifically,  we  show  that  there  exist  a  unique  fifth-order  RD  scheme  with  the 
generalized  trapezoidal  rule,  and  that  mathematically  becomes  sixth-order  RD  scheme.  Seeking  fifth-order 
accuracy  from  the  fourth-order  RD  scheme  given  in  Eq.  (58),  we  find  an  additional  constraint: 


CR  h 2 

+  =  (65) 

Then,  the  number  of  constraints  (four)  matches  the  number  of  unknown  coefficients  (four).  In  this  case, 
therefore,  there  exists  a  unique  solution: 


Gf  =  + 


lR 


CL  -  —I— h  R 
°2  "  +  60’ 


Cf  =  — 


C?  = 


60  ' 


(66) 


Interestingly,  the  above  unique  coefficients  satisfy  the  following  constraint  as  well,  which  is  the  requirement 
for  sixth-order  accuracy: 


cl 

4 


hR  +  C? 


h 2r 
30  ' 


(67) 
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Thus,  we  have  developed  a  sixth-order  RD  scheme  with  the  generalized  trapezoidal  rule.  We  remark  that 
these  coefficients,  unlike  the  ones  proposed  for  the  third-order  and  fourth-order  schemes,  are  unique.  With 
the  proposed  constants  for  the  sixth-order  RD  scheme,  the  source  term  evaluations  at  the  nodes  i  and  i  +  1 
become: 


Sl  =  Si  +  ^dxSi  +  ^dxxSi,  (68) 

5  oU 

SR  =  Si+1-f^fdxSi+1  +  f^dxxSi+1.  (69) 

Note  that  it  requires  only  the  first  and  second  derivatives  of  the  source  term,  and  no  higher-order  derivatives 
are  required. 

Similar  to  the  process  explained  for  the  fourth-order  scheme,  we  prove  the  order  of  accuracy  of  the  scheme 
using  the  above  proposed  source  discretization  by  first  evaluating  the  cell  residual  of  the  hyperbolic  system. 
Consider  the  cell  residual  of,  for  example,  the  first  equation  ($„ ): 

=  -a(ui+i-Ui)  +  iy(pi+1-pi)  +  ^Y(SL  +  SR), 

=  -a(ui+ 1  -  Ui )  +  v{pi+ 1  -  p^ 

+  _^_('S'i+i  +  Si)  — ^(ds-Si+i  —  dxSi)  +  Y^(dxocSi+ 1  +  dxxSi).  (70) 

Expanding  the  cell  residual  around  the  node  i.  we  obtain  the  truncation  error  of  the  first  equation  (after 
some  algebra)  as 


T.E.(dTUi)  = 

hn 

(- adxUi  +  vdxpi  +  Si)  +  —  ( -adxxUi  +  vdxxpi  +  dxSi) 

+ 

h2 

— —  (  Q'dxxx'U'i  V^xxxPi  H-  &xx&i) 

6 

+ 

h 3 

(  Cidxxxx'U'i  +  vdxxxxPi  H-  ^xxx^i) 

+ 

(  Q'^xxxxxW'i  H-  vdxxxxxPi  H”  ^xxxx^i) 

+ 

hR 

^2q  (  ^^xxxxxx^i  “1“  vdxxxxxxPi  ^xxxxx^i) 

21 

T  5040^  ^'^XXXXXXx'U'i  T  vdxxxXXXxPi  T  rjQ  ^XXXXXX-Si)  T  0(/l  ) 

=  0  +  0(h6).  (71) 

The  similar  result  is  obtained  for  the  second  equation,  and  therefore  is  not  repeated  here.  The  above 
truncation  error  analysis  reveals  that  this  powerful  sixth-order  RD  scheme  could,  in  practice,  produce  almost 

seventh-order  accurate  results,  if  the  sixth  derivative  of  the  source  term  becomes  small;  the  sixth-order  term  is 
only  due  to  the  presence  of  ( h^t/100800)dxxxxxxSi  in  the  last  term  of  Eq.  (71).  We  will  show  such  interesting 
results  in  the  following  section. 

We  emphasize  that  the  sixth-order  RD  scheme,  similar  to  the  third-order  RD  scheme,  requires  the 
evaluation  of  the  first  and  second  derivatives  of  the  source  terms.  However,  these  derivatives  are  now 
required  to  be  evaluated  with  higher-order  accuracy.  For  this  sixth-order  RD  scheme,  we  need  third-order 
and  second-order  accurate  evaluations  of  the  first  and  the  second  derivatives  of  the  source  terms,  respectively, 
which  can  be  obtained  by  a  cubic  fit.  This  makes  the  developed  sixth-order  scheme  slightly  more  expensive 
than  the  third-order  scheme.  Nevertheless,  the  proposed  sixth-order  scheme  is  exceptionally  simple  and 
affordable. 


VI.  Results 

In  this  section  we  present  the  results  in  three  categories:  1)  steady  advection-diffusion  equation  for  high 
Reynolds  (or  Peclet)  number  applications,  2)  unsteady  linear  advection-diffusion,  and  3)  unsteady  nonlinear 
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advection-diffusion  problems.  We  obtain  the  results  with  all  the  proposed  RD  schemes  and  compare  them 
with  the  second-order  RD  scheme.  The  order  of  accuracy  results  are  also  compared  and  presented  for  each 
example. 

For  all  cases,  we  employ  the  Newton-GS  solver  as  described  in  Section  III.  It  is  essentially  Newton’s 
method  for  the  second-order  scheme,  but  an  approximate  Newton  method  for  higher-order  schemes  because 
the  Jacobian  matrix  derived  from  the  second-order  scheme  is  used  for  all  higher-order  schemes.  For  unsteady 
problems,  the  same  solver  is  used  to  solve  the  implicit-residual  equations  or  equivalently  to  compute  the 
pseudo  steady  solution  at  every  physical  time  step. 

A.  Steady  Linear  Advection-Diffusion 

Consider  the  advection-diffusion  equation  in  a;  £  (0, 1)  with  u( 0)  =  0  and  m(1)  =  0: 

dtu  +  adxu  =  vdxxu  +  s(x),  (72) 

where 

s(x)  =  v-K2  s\h(ttx)  +  a7r  cos(7ra).  (73) 

The  above  problem  has  a  fixed  analytical  solution  of  sin(7ra;)  for  any  advection  and  diffusion  coefficients. 
We  solved  this  problem  using  the  hyperbolic  advection-diffusion  method  discussed  in  Section  II  with  the 
proposed  high-order  RD  schemes  (see  Fig.  5.  Note  that  we  intentionally  used  a  very  coarse  grid  to  show 
that  the  high-order  schemes  could  produce  a  much  more  accurate  solution  on  such  a  “bad”  grid.)  We 


(b)  p(x)  =  du/dx 


Figure  5.  Comparisons  between  the  exact  and  the  numerical  results  using  the  proposed  high-order  RD  schemes 
for  the  steady  linear  advection-diffusion  problem  {Re  =  1)  on  a  nonuniform  grid  with  N  =  10. 

used  ranges  of  nonuniform  grids  and  solved  the  problem  with  the  Newton-GS  method.  For  each  Newton 
iteration,  the  GS  relaxation  were  conducted  until  three  orders  of  magnitudes  reduction  is  achieved  for  the 
linear  system.  The  computations  were  continued  until  the  residuals  of  both  the  solution  u  and  the  solution 
gradient  p  were  reduced  by  eight  orders  of  magnitude.  The  convergence  results  are  shown  in  Table  1,  where 
the  convergence  of  the  GS  relaxation  is  clearly  0(N )  and  is  independent  of  the  order  of  accuracy  of  the  RD 
scheme.  Also,  the  solution  was  converged  with  the  same  number  of  GS  relaxations  and  Newton  iterations 
regardless  of  the  RD  scheme  order.  Note  that  the  solutions  were  converged  with  a  very  small  number  of 
Newton  iterations,  typically  less  than  ten;  this  is  exceptionally  remarkable  for  the  approximate  Jacobian 
(second-order)  formulation  for  higher-order  schemes. 

The  accuracy  of  the  proposed  RD  schemes  were  verified  by  computing  the  Li  =  Y^h=i(UiXact  ~  Ui)/N. 
These  results  are  shown  in  Fig.  6  for  the  third-,  fourth-,  and  sixth-order  hyperbolic  RD  schemes  (see  also 
Tables  2,  3,  and  4).  The  RD  schemes  developed  with  the  divergence  formulation  are  denoted  by  RD-D ,  while 
the  schemes  developed  with  the  generalized  trapezoidal  rule  are  denoted  with  RD-GT.  The  results  show  that 
all  the  RD  schemes  achieve  the  design  order  of  accuracy  for  both  the  solution  u  and  the  solution  gradient 
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Table  1.  Steady  linear  advection-diffusion  problem,  Re  =  1  (Convergence  criteria:  Residuals  <10  8.) 


Number  of  nodes 

RD  Scheme  Order 

GS  relaxations/Newton  iteration 

Newton  iteration 

High-Order  Technique 

High- Order  Technique 

RD-D 

RD-GT 

RD-D 

RD-GT 

3rd 

168 

169 

10 

10 

50 

4th 

168 

168 

10 

10 

6th 

- 

168 

- 

10 

3rd 

325 

325 

8 

8 

100 

4th 

325 

325 

8 

8 

6th 

- 

325 

- 

8 

3rd 

670 

670 

7 

7 

200 

4th 

670 

670 

7 

7 

6th 

- 

670 

- 

7 

3rd 

1015 

1015 

7 

7 

300 

4th 

1015 

1015 

7 

7 

6th 

- 

1015 

- 

7 

3rd 

1703 

1703 

7 

7 

500 

4th 

1703 

1703 

7 

7 

6th 

- 

1703 

- 

7 

6rd 

3416 

3416 

7 

7 

1000 

4th 

3416 

3416 

7 

7 

6th 

- 

3416 

- 

7 

Third-Order  RD  Schemes 


Fourth-Order  RD  Schemes 


Sixth-Order  RD  Scheme 


Figure  6.  L\  error  of  the  proposed  high-order  RD  hyperbolic  schemes  for  the  steady  linear  advection-diffusion 
problem. 


15  of  29 


American  Institute  of  Aeronautics  and  Astronautics 


Downloaded  by  NASA  LANGLEY  RESEARCH  CENTRE  on  August  11, 2014  I  http://arc.aiaa.org  I  DOI:  10.2514/6.2014-2090 


p.  It  can  be  seen  also  that  the  difference  between  the  two  versions  of  third-order  RD  schemes  is  minor,  and 
the  same  is  true  for  the  fourth-order  RD  schemes. 


Table  2.  Spatial  accuracy  with  the  third-order  RD-GT  scheme  for  the  steady  linear  advect ion-diffusion  problem 

(Re  =  1). 


Number  of  nodes 

Li  error  of  u 

Order 

L i  error  of  p 

Order 

10 

9.54E-03 

8.59E-02 

25 

4.63E-04 

3.30 

4.41E-03 

3.24 

50 

4.85E-05 

3.25 

4.63E-04 

3.25 

100 

5.46E-06 

3.15 

5.16E-05 

3.17 

200 

6.46E-07 

3.08 

6.04E-06 

3.09 

300 

1.88E-07 

3.04 

1.75E-06 

3.05 

500 

3.99E-08 

3.03 

3.71E-07 

3.04 

1000 

4.93E-09 

3.02 

4.57E-08 

3.02 

Table  3.  Spatial  accuracy  with  the  fourth-order  RD-D  scheme  for  the  steady  linear  advection-diffusion  problem, 
Re  =  1.  (Note:  fourth-order  RD-D  and  RD-GT  schemes  produce  almost  identical  result.) 


Number  of  nodes 

L\  error  of  u 

Order 

L i  error  of  p 

Order 

10 

8.38E-03 

6.36E-02 

25 

2.84E-04 

3.69 

1.82E-03 

3.88 

50 

1.74E-05 

4.03 

1.02E-04 

4.16 

100 

1.04E-06 

4.07 

5.77E-06 

4.15 

200 

6.23E-08 

4.06 

3.35E-07 

4.11 

300 

1.21E-08 

4.04 

6.44E-08 

4.07 

500 

1.55E-09 

4.02 

8.15E-09 

4.05 

1000 

9.30E-11 

4.06 

4.94E-10 

4.04 

B.  Steady  Boundary  Layer  Problem 

Consider  the  advection-diffusion  equation  in  x  £  (0, 1)  with  u( 0)  =  0  and  m(1)  =  1: 

dtu  + adxu  =  vdxxu  + s(x),  (74) 


where 


7T 

s(x)  =  —  [acos(7ra)  +  nu  sin(7r:r)] ,  Re  =  a/v. 


Re 


(75) 


This  is  a  boundary  layer  problem  with  a  non-trivial  steady  state  solution  in  the  diffusion  limit  as  a  result  of 
the  source  term  addition.  This  equation  develops  a  very  narrow  boundary  layer  near  the  right  boundary 
( x  =  1)  when  the  advection  term  becomes  dominant.  The  exact  steady  state  solution  to  this  problem  is 
given  by  (see  also  Fig.  7). 

p-Re  _  ( x-l)Re  i 

ue*act(x)  =  0-Re_l - +  777  sin(™)-  (76) 


Re 


We  chose  various  Re  values  ranging  from  1  to  106  and  solved  the  equation  on  nonuniform  grid  sizes 
up  to  100000  nodes.  Like  the  previous  problem,  the  solutions  were  obtained  with  the  Newton-GS  method. 
Within  each  Newton  iteration,  the  GS  relaxation  were  conducted  until  three  orders  of  magnitude  reduction  is 
achieved  for  the  linear  system,  and  the  residuals  of  both  the  solution  and  the  solution  gradient  were  reduced 
by  eight  orders  of  magnitude.  The  convergence  data  are  given  in  Table  5  for  the  proposed  high-order 
RD  schemes.  The  data  are  shown  for  the  high-order  schemes  developed  with  the  generalized  trapezoidal 
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Table  4.  Spatial  accuracy  with  the  sixth-order  RD-GT  scheme  for  the  steady  linear  advection-diffusion  problem 

(Re  =  l). 


Number  of  nodes 

Li  error  of  u 

Order 

L\  error  of  p 

Order 

10 

2.60E-03 

1.75E-02 

25 

3.34E-05 

4.75 

1.98E-04 

4.89 

50 

4.87E-07 

6.10 

2.87E-06 

6.11 

100 

4.32E-09 

6.82 

2.65E-08 

6.76 

200 

3.09E-11 

7.13 

1.78E-10 

7.22 

300 

5.80E-12 

4.13 

3.39E-11 

4.10 

500 

3.60E-12 

0.93 

1.97E-11 

1.06 

Re=10 


Re=10 


(a)  u(x) 


(b)  p(x)  =  du/dx 


Figure  7.  Comparisons  between  the  exact  and  the  numerical  results  using  the  proposed  high-order  RD  schemes 
for  the  steady  boundary  layer  problem  (Re  =  10)  on  a  nonuniform  grid  with  N  =  10. 
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rule.  Similar  results  were  obtained  with  the  RD  schemes  developed  with  the  divergence  formulation,  and 
therefore  not  shown.  But  in  general,  the  high-order  RD  schemes  with  the  generalized  trapezoidal  rule 
perform  better  than  the  ones  developed  with  the  divergence  formulation.  The  latter  third-order  RD-D 
scheme  not  only  produce  slightly  larger  errors  as  shown  in  Fig.  6,  but  also  encounter  some  convergence 
difficulties  (particularly  with  time-dependent  problems);  the  convergence  problem  seems  originated  from  the 
lack  of  diagonal  dominance  caused  by  the  particular  structure  of  the  source  term  discretization  as  shown 
in  Section  A.  We  will  therefore  mainly  focus  on  the  RD-GT  schemes  for  unsteady  cases.  The  0(N )  linear 


Table  5.  Steady  boundary  layer  problem  (Convergence  criteria:  Residuals  <  10  8.) 


logio  Re 

Number  of  nodes 

RD-GT  Scheme  Order 

GS  relaxations/Newton  iteration 

Newton  iteration 

3rd 

163 

8 

0 

50 

4th 

163 

8 

6th 

3rd 

324 

7 

0 

100 

4th 

324 

7 

6th 

3rd 

1647 

7 

0 

500 

4th 

1647 

7 

6th 

3rd 

178 

7 

1 

100 

4th 

178 

7 

6th 

3rd 

44 

7 

2 

100 

4th 

44 

7 

6th 

3rd 

42 

8 

3 

500 

4th 

43 

7 

6th 

3rd 

18 

9 

4 

1000 

4th 

19 

7 

6th 

3rd 

24 

9 

5 

10000 

4th 

21 

7 

6th 

3rd 

18 

9 

6 

100000 

4th 

19 

7 

6th 

dependency  of  the  GS  relaxations  on  the  grid  size  is  also  demonstrated,  which  is  a  consequence  of  solving  the 
advection-diffusion  equation  as  a  hyperbolic  system.  We  emphasize  that  this  is  remarkable  because  the  linear 
convergence  is  retained  for  any  irregular  grid  in  any  dimensions  (N  is  approximately  the  number  of  nodes 
in  each  coordinate  direction  in  two  and  three  dimensions)  as  demonstrated  in  Refs.  “  It  leads  to  orders 
of  magnitude  faster  convergence  in  comparison  with  conventional  methods  whose  convergence  is  typically 
0(N2)  as  discussed  also  in  the  previous  paper.  Furthermore,  the  proposed  high-order  RD  schemes  are 
extremely  efficient  as  the  solutions  were  obtained  with  only  a  small  number  of  Newton  iterations:  less  than 
10  iterations  to  reduce  the  residual  by  eight  orders  of  magnitude.  Moreover,  the  number  of  GS  relaxations 
and  Newton  iterations  are  essentially  independent  of  the  scheme  order.  Considering  the  fact  that  the  cost  of 
one  GS  relaxation  is  significantly  cheaper  than  one  Newton  iteration,  we  find  that  the  developed  high-order 
RD  schemes  are  extremely  powerful  and  efficient.  Finally,  as  in  the  previous  work,  we  remark  that  the 
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high-Re  cases  required  extremely  fine  grids  to  meet  the  well-known  requirement  on  the  mesh  Reynolds- 
number. J  If  desired,  the  computations  can  be  performed  on  substantially  coarser  grids  with  more  aggressive 
and  customized  grid  stretching.  However,  we  simply  refined  the  grid  to  meet  the  mesh  Reynolds-number 
requirement  because  our  method  is  powerful  enough  to  solve  the  problem  very  efficiently  (i.e. ,  less  than  10 
Newton  iterations)  even  on  such  dense  grids  for  all  third-,  fourth-,  and  sixth-order  schemes.  The  ability  to 
efficiently  solve  the  problem  on  highly  refined  grids  is  a  great  advantage  of  these  schemes. 

The  order  of  accuracy  of  the  proposed  RD  schemes  were  also  verified  for  this  problem.  Figures  8  shows 
the  L\  error  convergence  results,  where  h  is  the  representative  mesh  spacing  defined  by  h  =  /( N  —  1).  For 
discussion  purposes,  we  present  the  accuracy  plots  for  Re  =  1  and  Re  =  10;  similar  results  were  obtained  for 
other  Reynolds  numbers.  These  results  verify  the  order  of  accuracy  of  the  proposed  high-order  RD  schemes 
(i.e.  third-,  fourth-,  and  sixth-order)  for  all  the  variables  and  the  gradients  at  all  the  grid  nodes  including 
the  boundary  nodes  (see  also  Tables  6,  7,  and  8).  Note  that  when  the  differences  between  the  numerical 
results  and  exact  values  approach  the  machine  zero,  the  L\  slope  approaches  zeroth-order. 


Third-Order  RD  Schemes 

Fourth-Order  RD  Schemes 

Sixth-Order  RD  Schemes 

-l-U  (RD-GT) 

Hhu  (RD-GT) 

— | —  u  (RD-GT) 

4 

_0.p=ux  (RD-GT) 

4 

,0.p=ux  (RD-GT) 

4 

^.p=ux  (RD-GT) 

A  U  (RD-D) 

A  u  (RD-D) 

^  p=ux  (RD-D) 

p=ux  (RD-D) 

0 

0 

0 

_? 

, 

./-  /Or 

jf 

Slop.2  , 

jf  _4 

siop.2 

tf 

siop.2 

-6 

-6 

-6 

-B 

slop°3 

-8 

S1”P°3 

-8 

Slop.  3  s'  /jy 

-10 

-10 

„T. 

-10 

Slope  4  vsy 

-» 

Slope  bS*' T  t 

lof‘0h 

log'b 

lof‘0h 

(a)  third-order  (Re  =1) 

(b)  fourth-order  (Re  =  1) 

(c)  sixth-order  (Re  =  10) 

Figure  8.  L\  error  of  the  proposed  high-order  RD  hyperbolic  schemes  for  the  boundary  layer  problem  on 
nonuniform  grids. 


Table  6.  Spatial  accuracy  with  the  third-order  RD-GT  scheme  for  the  boundary  layer  problem  with  Re  =  1.0. 


Number  of  nodes 

Li  error  of  u 

Order 

Li  error  of  p 

Order 

10 

9.65E-03 

8.68E-02 

25 

4.70E-04 

3.30 

4.47E-03 

3.24 

50 

4.94E-05 

3.25 

4.70E-04 

3.25 

100 

5.56E-06 

3.15 

5.26E-05 

3.16 

200 

6.58E-07 

3.08 

6.16E-06 

3.09 

300 

1.91E-07 

3.06 

1.79E-06 

3.05 

500 

4.07E-08 

3.03 

3.79E-07 

3.04 

1000 

5.04E-09 

3.01 

4.66E-08 

3.02 

2000 

6.34E-10 

2.99 

5.79E-09 

3.01 

5000 

4.89E-11 

2.80 

3.68E-10 

3.00 

C.  Unsteady  Linear  Advection-Diffusion 

Consider  the  time-dependent  advection-diffusion  equation  in  a;  £  (0, 1) 


dtu  +  adxu  =  vdxxu.  (77) 

The  above  equation  with  the  following  initial  condition: 

u(x,t  =  0)  =  sin(«;x),  (78) 


19  of  29 


American  Institute  of  Aeronautics  and  Astronautics 


Downloaded  by  NASA  LANGLEY  RESEARCH  CENTRE  on  August  11, 2014  I  http://arc.aiaa.org  I  DOI:  10.2514/6.2014-2090 


Table  7.  Spatial  accuracy  with  the  fourth-order  RD-GT  schemes  for  the  boundary  layer  problem  with  Re  =  1.0 
Note  that  when  the  differences  between  the  numerical  results  and  exact  values  approach  the  machine  zero 
the  L\  slope  approaches  zeroth-order. 


Number  of  nodes 

Li  error  of  u 

Order 

L i  error  of  p 

Order 

10 

8.48E-03 

6.43E-02 

25 

2.88E-04 

3.69 

1.84E-03 

3.88 

50 

1.78E-05 

4.02 

1.04E-04 

4.14 

100 

1.06E-06 

4.07 

5.87E-06 

4.15 

200 

6.41E-08 

4.05 

3.42E-07 

4.10 

300 

1.25E-08 

4.04 

6.57E-08 

4.07 

500 

1.61E-09 

4.01 

4.35E-09 

4.05 

1000 

1.06E-10 

3.92 

4.98E-10 

4.06 

2000 

1.42E-11 

2.90 

1.42E-11 

3.40 

5000 

9.69E-12 

0.42 

4.75E-11 

0.01 

Table  8.  Spatial  accuracy  with  the  sixth-order  RD-GT  scheme  for  the  boundary  layer  problem  with  Re  =  10 


Number  of  nodes 

Li  error  of  u 

Order 

L\  error  of  p 

Order 

10 

4.81E-02 

5.54E-01 

25 

6.89E-05 

4.75 

6.19E-04 

4.89 

50 

4.54E-07 

6.11 

3.76E-06 

6.11 

100 

3.97E-09 

6.78 

3.21E-08 

6.73 

200 

3.56E-11 

7.59 

4.28E-10 

7.37 

300 

1.14E-11 

3.72 

1.40E-10 

4.23 

500 

9.23E-12 

0.05 

1.19E-10 

0.30 
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where  k  is  an  arbitrary  constant,  has  the  following  exact  solution  with  a  periodic  boundary  condition: 


uexact(x,  t )  =  e~K  vt  sin(n(x  —  at)). 

A  non-periodic  solution  also  exists  for  the  following  oscillatory  boundary  conditions 

u(0,t)  =  0, 

u(l,t)  =  U  cos(wf), 


(79) 


(80) 

(81) 


where  U  is  the  amplitude  of  the  oscillation  and  ui  is  the  frequency  of  the  oscillation  on  the  right  boundary. 
The  exact  solution  is  given  by 


uexact{x,t)  =  Real 


Al,2  — 


a  ±  \/ a2  +  4 iuv 
2v 


(82) 


where  i  =  \J — 1. 

We  solved  the  above  two  problems  with  the  first-order  hyperbolic  advection-diffusions  equation  given  as 
in  Eq.  (2)  and  (3).  For  each  physical  time,  we  reduced  the  residuals  by  two  orders  of  magnitude  before 
advancing  in  time.  During  each  time  step,  we  also  relaxed  the  linear  system  using  GS  relaxations  until  two 
orders  of  magnitude  reduction  in  the  linear  system  residuals  was  achieved  (see  Figs.  9  and  10.)  We  also  note 
that  more  residual  reduction  in  pseudo  time  may  be  necessary  on  more  complex  problems  but  two  orders  of 
magnitude  reduction  in  the  residuals  were  sufficient  for  the  problems  presented  here. 


(a)  u(x,t  =  1.0)  (b)  p(x,t  =  1.0)  =  du/dx 


Figure  9.  Time-dependent  linear  advection-diffusion  problem  (Re  =  33.33)  with  periodic  boundary  condition 
on  IV  =  10  uniform  nodes  (At  =  0.01.) 

We  examined  the  convergence  rate  of  these  problems  on  several  uniform  and  nonuniform  grid  systems. 
Given  in  Tables  9  and  10  are  the  average  numbers  of  GS  relaxations  per  Newton  iteration  obtained  over 
100  time  steps  for  the  periodic  and  the  oscillatory  boundary  condition  problems,  respectively.  Clearly  the 
convergence  rate  of  the  GS  relaxation  is  of  0(N ),  not  0(N2)  as  typical  for  numerical  methods  for  the 
advection-diffusion  equation.  Observe  that  for  most  grid  systems  only  two  Newton  iterations  were  sufficient 
to  obtain  accurate  solutions  regardless  of  the  scheme  order.  We  also  note  that  we  used  At  =  0.01  for  all 
the  grid  systems.  The  time  step  is  orders-of-magnitude  larger  than  that  required  for  conventional  explicit 
schemes,  which  is  limited  by  0(h2).  Of  course,  conventional  implicit  schemes  also  allow  unconditionally 
large  time  steps,  but  it  requires  0{N 2)  convergence  in  an  iterative  linear  solver  and  potentially  a  much 
larger  number  of  outer  iterations  as  well  if  the  exact  linearization  is  not  possible  and  Newton’s  method 
cannot  be  constructed.  Hence,  the  method  developed  here  has  two  major  advantages  over  conventional 
methods:  the  second  order  Jacobian  formulation  and  0(N)  iterative  convergence  in  the  linear  solver.  The 
latter  advantage  can  be  potentially  huge  with  increase  of  the  grid  system  as  the  speed-up  factor  is  O(N)  and 
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(a)  u(x,t  =  1.0)  (b)  p(x,t  =  1.0)  =  du/dx 


Figure  10.  Time-dependent  linear  advection-diffusion  problem  with  oscillatory  boundary  condition  (u  =  77t/2, 
ac  =  2-7T,  U  =  1,  v  =  1)  on  N  =  10  uniform  nodes  (At  =  0.01.) 


Table  9.  Unsteady  linear  advection-diffusion  problem  with  periodic  BC  (a  =  l,i/  =  0.03)  on  uniform  grids. 
Average  data  over  100  time  steps  are  given.  (Convergence  criteria:  Residuals  <  10— 2) 


Number  of  nodes 

RD-GT  Scheme  Order 

GS  relaxations/Newton  iteration 

Newton  iteration 

3rd 

5 

4 

25 

4th 

5 

4 

6th 

6 

5 

3rd 

24 

2 

100 

4th 

24 

2 

6th 

24 

2 

3rd 

33 

2 

300 

4th 

35 

2 

6th 

35 

2 

3rd 

55 

2 

500 

4th 

55 

2 

6th 

55 

2 

3rd 

116 

2 

1000 

4th 

116 

2 

6th 

116 

2 
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Table  10.  Unsteady  linear  advection-diffusion  problem  with  oscillatory  BC  (uj  =  77r/2,a  =  1.)  on  nonuniform 
grids.  Average  data  over  100  time  steps  are  given.  (Convergence  criteria:  Residuals  <  10~2) 


Number  of  nodes 

RD-GT  Scheme  Order 

GS  relaxations/Newton  iteration 

Newton  iteration 

3rd 

35 

4 

25 

4th 

37 

4 

6th 

45 

4 

3rd 

69 

4 

50 

4th 

72 

4 

6th 

72 

4 

3rd 

136 

4 

100 

4th 

142 

4 

6th 

142 

4 

3rd 

650 

4 

500 

4th 

680 

4 

6th 

680 

4 

3rd 

1283 

4 

1000 

4th 

1332 

4 

6th 

1332 

4 

thus  grows  for  finer  grids.  Note  also  that  the  second-order  Jacobian  formulation  is  the  advantage  of  the  RD 
scheme  over  finite  volume  schemes  where  the  compact  Jacobian  formulation  is  only  first-order.  The  results 
also  show  that  the  convergence  rate  is  the  same  for  all  the  developed  high-order  RD  schemes.  It  means  that 
the  only  cost  to  these  higher  order  schemes  is  the  evaluation  of  the  first  and  the  second  derivatives  of  the 
source  terms  and  in  the  case  of  the  fourth-order  scheme  with  the  generalized  trapezoidal  rule,  the  second 
derivatives  are  not  required. 

We  verified  the  order  of  accuracy  of  the  proposed  schemes  on  time-dependent  linear  problems  with 
consistent  space-time  discretization;  i.e.  third-order  with  BDF3,  fourth-order  with  BDF4,  and  sixth-order 
with  BDF6.  The  same  spatial  order  of  accuracy  can  be  observed  with  the  A-Stable  BDF2  with  small 
enough  time  steps  such  that  the  temporal  error  is  comparable  to  the  spatial  error  (i.e.  At2  ~  hm ,  where 
m  is  the  scheme  order) .  But  the  consistent  pair  of  BDF  and  spatial  discretization  allows  about  two  orders 
of  magnitude  reduction  in  the  number  of  time  steps  for  the  finest  grid  cases.  Note  that  time-accurate 
computations  are  started  by  BDF1  in  the  first  time  step  with  extremely  small  time  step  (e.g.  At  =  10-8), 
and  then  by  higher  order  BDFs  thereafter  with  much  larger  time  steps  (e.g.  At  =  0.01).  This  will  ensure 
the  order  of  accuracy  of  the  developed  higher  order  schemes  through  all  times.  We  remark  that  explicit  time 
stepping  is  not  available  for  time- accurate  computations  with  the  hyperbolic  system  method  (see  Ref.  for 
more  details.) 

Figure  11  shows  the  L\  error  convergence  for  the  above  time-dependent  linear  problem  with  the  oscillatory 
boundary  condition  (see  also  Tables  11,  12,  and  13),  where  clearly  the  order  of  accuracy  of  the  proposed 
schemes  are  verified  for  the  linear  time-dependent  problems.  The  results  were  obtained  at  t  =  1.0.  We  used 
the  same  At  among  the  third-,  fourth-  and  sixth-order  schemes  and  were  able  to  obtain  the  desired  order  of 
accuracy.  (Note  that  we  showed  the  RD  schemes  that  were  developed  with  the  generalized  trapezoidal  rule 
as  these  are  more  efficient  and  less  expensive  than  the  ones  developed  with  divergence  formulation  of  the 
source  terms.) 

D.  Unsteady  Nonlinear  Advection-Diffusion 

Consider  the  unsteady  nonlinear  viscous  Burgers  equation  with  an  unsteady  time-dependent  source  term: 

dtu  +  dxf  =  dx  (vux)  +  S(x,t),  a;  €(0,1),  (83) 
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Third-Order  RD  Scheme 


Fourth-Order  RD  Schemes 


Sixth-Order  RD  Scheme 


(a)  third-order  +  BDF3  (b)  fourth-order  +  BDF4  (c)  sixth-order  +  BDF6 

Figure  11.  Spatial  accuracy  of  the  proposed  high-order  RD-GT  hyperbolic  schemes  for  the  time-dependent 
linear  problem  with  oscillatory  BC  on  uniform  grids. 


Table  11.  Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC  (lj  =  7ir/2 ,a  =  1.)  using 
the  third-order  RD-GT  scheme  with  the  BDF3  time  discretization. 


Number  of  nodes 

At  (BDF3) 

L\  error  of  u 

Order 

L i  error  of  p 

Order 

10 

2.50E-03 

1.40E-04 

2.86E-04 

20 

2.50E-03 

9.06E-06 

3.95 

1.90E-05 

3.91 

50 

1.25E-03 

3.54E-07 

3.54 

7.38E-07 

3.54 

100 

5.00E-04 

2.57E-08 

3.78 

4.71E-08 

3.97 

200 

2.50E-04 

2.59E-09 

3.31 

4.40E-09 

3.42 

Table  12.  Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC  (cj  =  7tt/2 ,a  =  1.)  using 
the  fourth-order  RD-GT  scheme  with  the  BDF4  time  discretization. 


Number  of  nodes 

At  (BDF4) 

L\  error  of  u 

Order 

L i  error  of  p 

Order 

10 

2.50E-03 

1.32E-04 

2.83E-04 

20 

2.50E-03 

7.22E-06 

4.19 

1.63E-05 

4.12 

50 

1.25E-03 

1.70E-07 

4.09 

4.01E-07 

4.04 

100 

5.00E-04 

1.04E-08 

4.03 

2.51E-08 

4.00 

200 

2.50E-04 

6.78E-10 

3.94 

1.64E-09 

3.94 

Table  13.  Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC  (uj  =  77r/2,a  =  1.)  using 
the  sixth-order  RD-GT  scheme  with  the  BDF6  time  discretization. 


Number  of  nodes 

At  (BDF6) 

L\  error  of  u 

Order 

L i  error  of  p 

Order 

10 

2.50E-03 

3.77E-06 

2.29E-05 

20 

2.50E-03 

7.19E-08 

5.71 

3.27E-07 

6.13 

50 

1.25E-03 

5.81E-09 

6.20 

3.70E-08 

5.37 

100 

5.00E-04 

3.62E-10 

5.43 

1.34E-09 

6.50 
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where  /  =  it2/ 2,  v  =  u,  and 


(84) 


S(x,t)  =  dtue  +  -dx(ue)2  -  dx{uepe ), 

where  pe  =  dxue.  The  source  term  has  been  generated  by  the  following  function: 


u'(x,t)  =  Real  ( sillh(l/gZ;)ce‘",)  +C,  C  >  1, 

\  sinh(-\/iw/i/)  / 

so  that  it  is  the  exact  solution  to  Eq.  (83)  with  the  boundary  conditions  defined  as 


u(0,t)  =  C , 

u(l,t)  =  C  +  (7cos(wf), 


(85) 


(86) 

(87) 


where  ui  is  the  frequency  of  the  oscillation  on  the  right  boundary,  and  U  is  the  amplitude  of  the  oscillation. 
We  note  that  the  constant  C  must  be  greater  than  1  in  order  for  the  diffusion  coefficient  to  be  positive.  We 
solved  this  time-dependent  nonlinear  advection-diffusion  equation  with  the  following  equivalent  first-order 
hyperbolic  system  (see  Ref.  for  more  details): 


dTu  +  dx  (w2)  =  dxp  —  dtu  +  S(x,t), 

(88) 

T 

—  dTp  =  (dxu-p/v). 

(89) 

For  this  nonlinear  unsteady  problem,  the  manufactured  source  term  contains  terms  that  are  already  in  the 
divergence  form;  i.e.  the  residual  evaluation  of  these  terms  are  exact.  The  dtue  term  is  the  only  term  in  the 
manufactured  source  term  with  a  non-exact  residual  evaluation.  The  BDF  discretization  of  the  dtu  term 
in  Eq.  (88)  will  not  be  in  the  divergence  form  and  therefore  will  not  have  an  exact  residual  evaluation.  In 
addition,  the  second  equation  has  a  nonlinear  source  term,  p/ v  (note  that  here  v  =  u).  We  obtained  the 
high-order  results  by  evaluating  all  of  these  non-exact  residuals  using  the  proposed  techniques.  Newton 
iterations  are  taken  to  be  converged  when  the  overall  residuals  are  dropped  by  eight  orders  of  magnitude. 
For  each  Newton  iteration,  we  relaxed  the  linear  system  until  the  residuals  are  reduced  by  two  orders  of 
magnitude. 

The  O(N)  convergence  rate  of  the  GS  relaxations  was  once  again  achieved  for  the  time-dependent  nonlin¬ 
ear  hyperbolic  advection-diffusion  system.  This  is  given  in  Table  14,  where  the  average  number  of  iterations 
were  obtained  over  1000  time  steps  (over  17  periods).  Note  also  that  the  high-order  RD  schemes  converged 
with  only  ten  Newton  iterations  with  the  compact  second-order  Jacobian  formulation.  It  is  remarkable  that 
such  a  rapid  convergence  is  achieved  for  all  high-order  schemes  with  the  second-order  Jacobian. 


Table  14.  Unsteady  nonlinear  advection-diffusion  problem  (u  =  a  =  u)  with  oscillatory  BC  (U  =  1,  C  =  2, 
lo  =  7tt/2)  on  nonuniform  grids.  Average  data  over  1000  time  steps  are  given.  (Convergence  criteria:  GS 
Relaxation  <  1CU2;  Newton  residuals  <  10~s). 


Number  of  nodes 

RD-GT  Scheme  Order 

GS  relaxations/Newton  iteration 

Newton  iteration 

3rd 

435 

10 

50 

4th 

430 

10 

6th 

431 

10 

3rd 

879 

10 

100 

4th 

868 

10 

6th 

864 

10 

3rd 

1772 

10 

200 

4th 

1749 

10 

6th 

1737 

10 

Shown  in  Fig.  12  are  the  order  of  accuracy  plots  obtained  for  this  unsteady  advection-diffusion  problem 
on  series  of  nonuniform  grids  (see  also  Tables  15,  16,  and  17).  The  results  confirm  the  high-order  accuracy 
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of  developed  RD  schemes  for  unsteady  nonlinear  problems  on  nonuniform  grids.  We  also  observe  that  our 
proposed  sixth-order  RD  scheme,  as  discussed  in  the  previous  section,  have  in  fact  produced  seventh-order 
accurate  results  demonstrating  its  power  characteristics  in  efficiently  providing  very  high  accurate  solutions 
and  gradients  (see  Table  17  for  more  details).  We  note,  however,  that  the  seventh-order  accuracy  is  not  a 
general  feature  of  the  sixth-order  scheme.  Although  it  is  not  clear  from  the  exact  solution,  it  could  be  due 
to  the  sixth-order  derivative  of  the  source  term  in  the  leading  truncation  error  being  vanishingly  small. 


Third-Order  RD  Scheme 


Fourth-Order  RD  Schemes 


Sixth-Order  RD  Scheme 


(a)  third-order  +  BDF3 


-e-p(=ux> 


loS.oh 

(b)  fourth-order  +  BDF4 


(c)  sixth-order  +  BDF6 


Figure  12.  Spatial  accuracy  of  the  proposed  high-order  RD-GT  hyperbolic  schemes  for  the  time-dependent 
nonlinear  problem  (y  =  a  =  u)  with  oscillatory  BC  ( U  =  1,  C  =  2,  oj  =  77r/2)  on  nonuniform  grids. 


Table  15.  Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (y  —  a  =  u)  with  oscillatory  BC  ( U  =  1, 
C  =  2,  uj  =  7tt/2)  using  the  third-order  RD-GT  scheme  and  BDF3  time  discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF3) 

L\  error  of  u 

Order 

L\  error  of  p 

Order 

25 

1.25E-03 

6.83E-05 

3.59E-04 

50 

1.25E-03 

3.68E-06 

4.21 

1.79E-05 

4.33 

100 

1.00E-03 

2.33E-07 

3.98 

1.12E-06 

4.00 

200 

5.00E-04 

1.70E-08 

3.78 

8.23E-08 

3.77 

500 

2.50E-04 

9.58E-10 

3.14 

5.03E-09 

3.05 

Table  16.  Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (y  =  a  =  u)  with  oscillatory  BC  (U  =  1, 
C  =  2,  to  =  77t/2)  using  the  fourth-order  RD  schemes  and  BDF4  time  discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF4) 

L\  error  of  u 

Order 

L\  error  of  p 

Order 

25 

1.25E-03 

6.75E-05 

3.49E-04 

50 

1.25E-03 

6.21E-07 

4.27 

2.71E-06 

4.23 

100 

1.00E-03 

1.85E-07 

4.24 

7.68E-07 

4.40 

200 

1.00E-03 

1.07E-08 

4.11 

3.69E-08 

4.38 

500 

5.00E-04 

3.14E-10 

3.85 

6.33E-10 

4.44 
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Table  IT.  Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (v  =  a  =  u)  with  oscillatory  BC  ( U  =  1, 
C  =  2,  uj  =  77t/2)  using  the  sixth-order  RD-GT  scheme  and  BDF6  time  discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF6) 

L\  error  of  u 

Order 

L\  error  of  p 

Order 

50 

1.25E-03 

7.99E-08 

5.61E-08 

75 

5.00E-04 

5.42E-09 

6.60 

3.69E-09 

6.71 

100 

5.00E-04 

7.54E-10 

6.85 

5.09E-09 

6.91 

200 

5.00E-04 

4.30E-11 

7.06 

2.83E-09 

7.11 

VII.  Conclusions 

In  this  paper,  we  have  developed  a  series  of  efficient  high-order  Residual-Distribution  (RD)  schemes 
for  general  advection-diffusion  problems.  Third-  and  fourth-order  RD  schemes  were  developed  with  the 
divergence  formulation  of  the  source  term.  These  schemes  are  very  economical  as  they  require  only  the 
evaluation  of  the  first  and  second  derivatives  of  the  source  term,  and  not  of  the  solution.  The  first  and 
second  derivatives  need  to  be  second  and  first  order  accurate,  respectively,  on  arbitrary  grids,  and  they 
are  obtained  by  a  compact  quadratic  fit.  Third-,  fourth-,  and  sixth-order  RD  schemes  are  also  developed 
with  a  generalized  trapezoidal  rule.  The  third-order  schemes  require  the  evaluation  of  the  first  and  second 
derivatives  of  the  source  term,  which  must  be  second  and  first  order  accurate,  respectively.  On  the  other 
hand,  the  fourth-order  scheme  requires  only  the  second-order  accurate  gradients,  and  therefore  is  the  least 
expensive  scheme  among  all  the  developed  high-order  schemes  in  this  paper.  All  third-  and  fourth-order 
schemes  are  constructed  within  a  five-point  stencil  in  the  interior,  a  four-point  stencil  at  the  nodes  adjacent 
to  the  boundary,  and  a  three-point  stencil  at  the  boundary  nodes.  For  the  sixth-order  scheme,  the  evaluation 
of  the  first  and  second  derivatives  of  the  source  term  is  required,  and  they  must  be  third-  and  second-order 
accurate,  which  is  achieved  by  a  cubic  fit.  It  results  in  the  stencil  extension  by  one  or  two  nodes  in  the  nodal 
residual.  In  addition,  the  analysis  of  the  proposed  sixth-order  RD  scheme  as  well  as  the  results  presented 
here  show  that  this  high-order  RD  scheme  can,  in  some  cases,  produce  seventh-order  results  on  nonuniform 
grids.  An  implicit  steady  solver  is  constructed  based  on  the  Jacobian  derived  from  the  compact  second-order 
RD  scheme.  It  is  demonstrated  that  the  solver  achieves  a  rapid  convergence  like  Newton’s  method  for  all 
high-order  schemes  despite  the  fact  that  the  Jacobian  is  not  exact.  Specifically,  it  requires  only  a  small 
number  of  Newton  iterations,  typically  less  than  10,  for  both  steady  and  unsteady  problems,  even  for  highly 
refined  grids,  up  to  100000  nodes.  We  have  demonstrated  also  that  all  of  these  high-order  schemes  are 
capable  of  producing  both  high-order  accurate  solution  and  gradient  on  nonuniform  grids  very  efficiently  by 
less  than  ten  Newton  iterations. 

The  study  presented  in  this  paper  should  be  of  interest  to  researchers  working  on  finite-volume  schemes 
because  the  RD  scheme  is  known  to  be  equivalent  to  the  finite-volume  scheme  in  one  dimension  (see  for 
example  Ref.  ).  Specifically,  a  second-order  upwind  RD  scheme  is  equivalent  to  the  first-order  finite- volume 
scheme  with  a  special  form  of  source  term  discretization.  '1  It  implies  that  the  developed  high-order  RD 
schemes  may  also  be  implemented  in  the  form  of  first-order  finite-volume  schemes  with  special  source  term 
discretization  formulas.  The  resulting  finite-volume  schemes  will  be  different  from  many  other  finite-volume 
schemes  in  that  they  do  not  require  computations  of  solution  gradients. 

The  developed  high-order  schemes  could  well  bring  significant  improvements  to  the  numerical  methods 
for  practical  problems  such  as  material  thermal  response  calculations  of  thermal  protection  systems  of  at¬ 
mospheric  entry  vehicles,'  “  and  the  experimental  aeroheating  data  reduction,'  ’  which  are  based  on 
one-dimensional  analyses  and  routinely  used  in  industries  (e.g.  NASA).  A  particularly  useful  scheme  would 
be  the  fourth-order  scheme  based  on  the  generalized  trapezoidal  rule  (RD-GT)  because  it  requires  only 
the  second-order  accurate  gradients  of  the  source  term.  Application  to  these  practical  problems  should  be 
undertaken  and  is  left  as  future  work. 

Extensions  to  higher  dimensions  are  highly  desired.  To  extend  the  developed  high-order  schemes  to  higher 
dimensions,  it  is  necessary  to  develop  a  high-order  quadrature  formula  for  integrating  the  flux  divergence 
term,  which  is  exact  in  one  dimension  but  not  in  higher  dimensions.  For  the  source  term  discretization,  the 
divergence  formulation  can  be  extended  relatively  straightforwardly  while  a  discretization  formula  such  as  the 
generalized  trapezoidal  rule  remains  to  be  found.  We  also  remark  that  the  BDF2  is  A-stable  but  higher-order 
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BDFs  are  not  and  may  encounter  instability  in  practical  problems  for  complex  equations  and  geometries  in 
higher  dimensions.  Alternative  implicit  time  integration  methods,  such  as  space-time  methods,  may  turn 
out  to  be  useful  or  necessary  for  such  applications.  Finally,  extensions  to  more  complex  nonlinear  equations 
such  as  the  Navier-Stokes  equations  remain  as  a  challenge.  For  the  compressible  Navier-Stokes  equations, 
the  complete  eigen-structure  of  the  whole  system  has  yet  to  be  found.  The  construction  of  the  upwind  RD 
scheme  based  on  a  single  hyperbolic  system,  as  presented  in  this  paper,  is  a  challenge.  To  overcome  the 
difficulty,  a  simplified  approach  has  been  proposed  and  demonstrated  for  a  finite- volume  scheme  in  Ref., ' 
which  is  based  on  the  independent  treatment  of  the  inviscid  and  viscous  terms.  A  similar  approach  may 
become  necessary  for  the  RD  schemes. 
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